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ABSTRACT 

We present a hadronic model for gamma-ray production in the microquasar 
LS I +61 303. The system is formed by a neutron star that accretes matter from 
the dense and slow equatorial wind of the Be primary star. We calculate the 
gamma-ray emission originated in pp interactions between relativistic protons in 
the jet and cold protons from the wind. After taking into account opacity effects 
on the gamma-rays introduced by the different photons fields, we present high- 
energy spectral predictions that can be tested with the new generation Cherenkov 
telescope MAGIC. 

Subject headings: gamma-rays: observations — gamma-rays: theory — X-ray: 
binaries — stars: individual (LS I +61 303) 
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1. Introduction 

LS I +61 303 is a Be/X-ray binary that presents unusually strong and variable radio 
emission (Gregory & Taylor 1978). The X-ray emission is weaker than in other objects of the 
same class (e.g. Greiner & Rau 2001) and shows a modulation with the radio period (Paredes 
et al. 1997). The most recent determination of the orbital parameters (Casares et al. 2005) 
indicates that the eccentricity of the system is 0.72 ± 0.15 and that the orbital inclination is 
~ 30°±20°. The best determination of the orbital period (P = 26.4960±0.0028) comes from 
radio data (Gregory 2002). The primary star is a B0 V with a dense equatorial wind. Its 
distance is ~ 2 kpc. The X-ray/radio ourtbursts are triggered 2.5-4 days after the periastron 
passage of the compact object, usually thought to be a neutron star. These outbursts can 
last until well beyond the apastron passage. 

Recently, Massi et al. (2001) have detected the existence of relativistic radio jets in 
LS I +61303, which makes of it a member of the microquasar class. Microquasars are 
thought to be potential gamma-ray sources (Paredes et al. 2000, Kaufman Bernado et 
al. 2002, Bosch-Ramon et al. 2005a) and, in fact, LS I +61303 has long been associated 
with a gamma-ray source. First with the COS-B source CG135+01, and later on with 3EG 
J0241+6103 (Gregory & Taylor 1978, Kniffen et al 1997). The gamma-ray emission is clearly 
variable (Tavani et al. 1998) and has been recently shown that the peak of the gamma-ray 
lightcurve is consistent with the periastron passage (Massi 2004), contrary to what happens 
with the radio/X-ray emission, which peaks after the passage. 

The matter content of microquasars jets is unknown, although in the case of SS 433 
iron X-ray line observations have proved the presence of ions in the jets (Kotani et al. 1994, 
1996; Migliari et al. 2002). In the present paper we will assume that relativistic protons 
are part of the content of the observed jets in LS I +61 303 and we will develop a simple 
model for the high-energy gamma-ray production in this system, with specific predictions for 
Cherenkov telescopes like MAGIC. We emphasize that our model is not opposed, but rather 
complementary to pure leptonic models as those presented by Bosch-Ramon & Paredes 
(2004) and Bosch- Ramon et al. (2005a), since the leptonic contribution might dominate 
at lower gamma-ray energies and after the periastron passage. In the next section we will 
describe the basic features of the model, and then we will present the calculations and results. 

2. General picture 

A hadronic model for the gamma-ray emission in microquasars with early-type compan- 
ions has been already developed by Romero et al. (2003). This model, however, is limited to 
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the simple case of a massive star with a spherically symmetric wind and a compact object in 
a circular orbit. Here we will consider a B-type primary with a wind that forms a circumstel- 
lar outflowing disk of density p w (r) = p (r/R*)~ n (Gregory & Neish, 2002). The continuity 
equation implies a wind velocity of the type v w = v (r / R*) n ~ 2 . We will consider that the 
wind remains mainly near to the equatorial plane, confined in a disk with half-opening angle 
= 15°, with n = 3.2, p = 10 -11 g cm" 3 , and v = 5 km s" 1 (Marti and Paredes 1995). 

The modeled properties of the system will be expressed in terms of the orbital phase 
ip (ip — 0.23 at the periastron passage according to the latest determination by Casares et 
al. 2005) which is related to the separation between the stars by r(ip) = a(l — e 2 )/[l — 
e cos(27r(?/> + 0.73))], where a is the semi-major axis of the orbit and e the eccentricity. The 
wind accretion rate onto the compact object of mass M c can be estimated as: 

it^y, (i) 

^rel 

where t> re i is the relative velocity between the neutron star (moving in a Keplerian orbit) and 
the circumstellar wind, assumed to be flowing radially on the equatorial plane. 

Following the basic assumption of the jet-disk symbiosis model (Falcke & Biermann 
1995) we will assume that the accretion rate is coupled to the kinetic jet power by: 

Qi = %M c c 2 , (2) 

where q$ ~ 0.1 is the coupling constant. Most of this power will consist of cool protons that 
are ejected with a macroscopic Lorentz factor T ~ 1.25 (Massi et al. 2001). Only a small 
fraction q- el ~ 10~ 3 is in the form of relativistic hadrons. The relativistic jet is confined by 
the pressure of the cold particles (P co id > -Prei), which expand laterally at the local sound 
speed. 

The jet axis, z, will be assumed normal to the orbital plane. The jet will be conical, 
with a radius Rj(z) = z(Ro/z ), where z is the injection point and R is the initial radius 
of the jet. We will adopt Zq = 10 7 cm and R$ = zq/10 as reasonable values (see Romero 
et al. 2003 and Bosch-Ramon et al. 2005a, who deals with similar jets for additional 
details). The relativistic proton spectrum will be a power law N'(E') = K p E'~ a , valid for 
E' p mm < E' < E' p max (in the jet frame). The corresponding relativistic proton flux will be 
Jp(E') = (c/ ' Atx)N' p {E' p ). Since the jet expands in a conical way, the proton flux evolves with 
z as: 2 

m) = (f ) E>r, (3) 

where it is implicit the assumption of the conservation of the number of particles (see Ghis- 
ellini et al. 1985), and a prime refers to the jet frame. Using relativistic invariants, it can 
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be shown that the proton flux, in the lab (observer) frame, becomes (e.g. Purmohammad & 
Samimi 2001): 



J P (E P , 9) 



cK fz Y r-«+! (E p - foyJEl -m 2 c 4 cos 9)' 



An V z J 



sin 2 9 + P cos 9 - 



1/2- 



(4) 



In this expression, T is the jet Lorentz factor, 9 is the angle subtended by the proton velocity 
direction (which will be roughly the same as that of the emerging photon) and the jet 
axis (notice that then 9 ~ 9 ^, s ), and (3^ is the bulk velocity in units of c. We will make 
all calculations in the lab frame, where the cross sections for pp interactions have suitable 
parametrizations. 

The number density no' of particles flowing in the jet at Rq, and the normalization 
constant K can be determined as in Romero et al. (2003). In the numerical calculations 
of the next section we have considered E'™ ax = 100 TeV, E'™ m = 1 GeV, T = 1.25, and, 
a = 2.2 (see the list of the assumed parameters in Table 1). The assumed maximum energy 
is consistent with the jet size and shock acceleration with an efficiency ~ 0.01 — 0.1. 

The matter from the wind can penetrate the jet from the side, diffusing into it as long 
as the particle gyro-radius is smaller than the radius of the jet. This imposes a constraint 
onto the value of the magnetic field in the jet: B jCt > E k /(eR ), where E k = m p t> 2 cl /2. For 
the periastron passage (E k maximum) results Bj Ct > 2.8 10~ 6 G, which is surely satisfied. 
However, some effects, like shock formation on the boundary layers, could prevent some 
particles from entering into the jet. Given our ignorance of the microphysics involved, we 
adopt a parameter f p that takes into account particle rejection from the boundary in a 
phenomenological way. In a conservative approach, we will adopt f p ~ 0.1 . 

Some of the particles entering the jet flow would be immediately accelerated to the jet 
velocity (by Coloumb interactions or wave-particle interactions). As a consequence, the jet 
should be slowed down during its motion through the equatorial wind. However, it is a fact 
that the jet survives this interaction since it is seen at radio wavelengths far beyond the 
wind region, up to distances of ~ 400 AU (Massi et al. 2001, 2004). Since the bulk velocity 
seems not to be very high (Massi et al. 2001) and hence its change does not affect seriously 
the calculations of the gamma-ray emissivity, we will neglect, in what follows, the effects of 
a macroscopic deceleration. The reader interested in the case of the hadronic gamma-ray 
emission of a jet slowed down to rest by the effects of the wind and the resulting standing 
shock wave as the major source of radiation is referred to the recent treatment presented by 
Romero & Orellana (2005). 
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In Figure 1 we show a sketch of the general situation and in Figure 2 we show the orbit 
of the system and the corresponding phases. 

3. Gamma-ray emission 

Relativistic protons in the jet will interact with target protons in the wind through 
the reaction channel p+p^p+p+ ^on + £„.±(7r + + 7r~), where £ n is the corresponding 
multiplicity. Then pion decay chains will lead to gamma-ray and neutrino emission. The 
differential gamma-ray emissivity from 7r°-decays can be expressed as (e.g. Aharonian & 
Atoyan 1996): 

g 7 (£ 7 , 9) = 47ra pp (£ p )— ^ J p (£ 7 , 9) rj A: (5) 

where Z^l^ is the so-called spectrum-weighted moment of the inclusive cross-section (see, 
for instance, Gaisser 1990). J P {E^) is the proton flux distribution (4) evaluated at E — E 1 . 
The cross section a pp (E p ) for inelastic p — p interactions at energy E p m Q^oE 7 /K, where 
K ~ 0.5 is the inelasticity coefficient and = l.l^p/GeV) 1 / 4 , can be represented for 
E p > 1 GeV by 

a pp (E p ) ss 30 x [0.95 + 0.06 log (E p /GeV)\ (mb). 

Finally, the parameter tja takes into account the contribution from different nuclei in the wind 
and in the jet (for standard composition of cosmic rays and interstellar medium t/a ~ 1.4). 

The spectral energy distribution is: 

L 7 (£ 7 , 0)=E% [ ra(r') g 7 (£ 7 , 9) d 3 f>, (6) 
Jv 

where V is the interaction volume between the jet and the circumstellar disk. The particle 
density of the wind that penetrates the jet is n(r) ~ f p p w (r)/m p . 

In our calculations, we adopt a viewing angle of 9 = 30° in accordance with the average 
value given by Casares et al. (2005). In Figure 3 we show a 3-D plot that shows the evolution 
of the gamma-ray spectral energy distribution as a function of the orbital phase. Other 
two plots in this figure show cuts at both the periastron and apastron, and the luminosity 
evolution with the orbital phase at 100 GeV. In both cases we show the unabsorbed (dashed 
lines) and the absorbed (continuum lines) curves. This absorption is discussed in the next 
section. 

At the periastron passage the unattenuated luminosity is ~ 10 33 erg s -1 . We can make 
a simple order-of-magnitude estimate of this value. The accretion rate at the periastron 
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is ~ 3 x 10 17 g s _1 . This means that the total power in relativistic protons should be 
Qj el = 10~ 3 M c c 2 ~ 2.8 x 10 35 erg s" 1 . The density of the stellar wind at the injection point 
of the jet is n ~ 4 x 10 11 cm~ 2 and the cross section for protons of E p ~ 1 TeV, a pp ~ 34 mb. 
Hence, the mean free path of the protons results \ pp ~ 8.3 x 10 13 cm. The thickness of the 
region of the disk traversed by the jet is Az ~ r perias tan 15° ~ 4.4 x 10 11 cm. Consequently, 
we can approximate the gamma-ray luminosity by: 

L 7 = 2UQ? (1 - e ~ Az ^) , (7) 

where f n ~ 0.2 is the fraction of the energy of the leading proton that goes into neutral pions 
and hence into gamma-rays. With a simple substitution into Eq. (7) we get L 7 ~ 6.6 x 10 32 
erg s -1 , in good agreement with the detailed numerical calculations presented in Fig. 3. 



4. Opacity 

The optical depth for a photon with energy E 7 , which in this case depends upon the 
direction observed, can be estimated as 



poo poo 

t(p,Et)= / n v ^(E vh ,p')cr e - e +(E vh ,E 1 )dp'dE vhl (8) 

^E min (E 7 ) Jp 

where E p ^ is the energy of the ambient photons, n p h(-E p h,p) is their density at a distance p 
from the neutron star, and cr e - e +(E ph , E 7 ) is the photon-photon pair creation cross section 
given by: 



2^ 2 ~ 2) + (3 - £ 4 ) In 



1+j 
1-C 



(9) 



where r is the classical radius of the electron and 



(m e c 2 ) 2 
EphEj 



1 - 



1/2 



(10) 



In Eq. (8), E m [ n is the threshold energy for pair creation in the ambient photon field. This 
field can be considered as formed by two components, one from the Be star and the other 
from the hot accreting matter impacting onto the neutron star: n p h = n ph l + n p h,2- Here, 



n p h,i(£ P h,p) 



Rl 



7rB(E ph ) \ 

lie Eph J p 2 + r 2 — 2pr sin 9 



(11) 
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is the black body emission from the star, with 

B{Eph) ~ (hey (e^/^« - 1) (12) 

and T cff = 22500 K (Marti & Paredes 1995). The separation r between the stars is again 
variable with the phase angle ip. 

The emission from the heated matter can be approximated by a Bremsstrahlung spec- 
trum: 

L E~^ 

w^ph, p) = 47rcp2 ^ ph /L- off for E ^ * 1 keV > ( 13 ) 

where Lx is the total luminosity in hard X-rays and £ , cut - fr ~ 100 keV. The photon index 
of the hard X-rays is taken to be within the range published by Greiner & Rau (2001), which 
was observationally determined. Lx is also constrained by observations, being Lx ~ 10 34 
erg s _1 (Paredes et al. 1997). Notice that no bump due to a putative accretion disk has 
been observed at X-rays, so we neglect this contribution. 

As an example, Figure 4 shows the dependence of the optical depth r with the energy 
of the 7-rays and its variation along the z axis for the observer at # bs = 30°. From detailed 
versions of this plot, we find that for photons of E 1 = 100 GeV significant absorption occurs 
mostly between ip — 0.1 and ip = 0.5. The optical depth remains well below the unity along 
the whole orbit for photons of energies E 1 < 30 GeV and E 1 > 2 TeV. 



5. Secondary electron-positron pairs and synchrotron emission 

Secondary pairs are produced by the decays of charged pions and muons, as well as by 
photon-photon interactions. The main reactions that lead to charged pions are: 

P + P -> p + p + ^07l° + ^±(7T + + 7T-) (14) 

p + p -> p + n + n + + X (15) 
p + p -> 2n + 2?r + + X (16) 

where n is a neutron, X stands for anything (neutral) else, and the charged pion multiplicity 
is £„.± PS 2(i£p/GeV) 1 / 4 . The neutrons have a proper lifetime of 886 ± 1 s and since they 
move at ultrarelativistic speed can escape from the source, decaying at considerable distances 
(Eichler & Wiita 1978). On the contrary, pions decay into the jet trough 7r ± — ^ /x ± H- and 
jj^ — > e ± + v + 77. For an injection proton spectrum given by Eq. (3) with a = 2.2, we have 
that the pion spectrum (in the jet's system) will be a power-law J' n± (E' w± ) = K v ± E'~± n , with 
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w 



a n ~ 2.3. The electron-positron distribution mimics this power law (Ginzburg & Syrovatskii 
1964, Dermer 1986): 

J> e± {E' e± )=K^ e± E>- e «\ (17) 

with 

A „ e± = f^r -1 2 <°* + 5 > - ^, us) 

Vm e / a±(a± + 2)(a± + 3) 

and 0:^ = 0;^. 

The energy density of pion-generated pairs along the jet at the periastron passage can 
be calculated as: 

We± = J(4rr/c)E' e ±f e ±(E' e ±)dE' e ±, (19) 

where J' e ±(E' e ±) takes into account all the contributions from z to -2 max - Integrating we get 
w^ e ± ~ 3 x 10 9 erg cm' 3 . 

We can compare the energy density of pairs from the charged pion decays with that of 
the pairs produced by direct gamma-ray absorption. The total luminosity of these pairs is: 

L e± = L;(l-e- T ). (20) 

Then, using the opacity calculated in the previous section, the pair energy density results 

w^ e± ~ (21) 

At the periastron passage, we get w 11 ^ e ± ~ 3.7 x 10 9 erg cm" 3 . Hence, the pair injection 
from the photon-photon annihilation is similar to that of pion decay. In what follows we will 
evaluate the spectrum of these particles using the approximation derived by Aharonian et 
al. (1983), which is in excellent agreement with the more detailed calculations (exact to 2nd 
order QED) presented by Bottcher & Schlickeiser (1997). 

The differential pair injection rate is given by (Bottcher & Schlickeiser 1997): 



3 [°° iV 7 (6 7 ) f°° n ph (e ph ) 
n e±(l) = ^ca T / de 7 — / de ph x 

6Z J-i e % r C n 



e 7 

.2 



7 J 47(e7 - 7 ) Ph 



7(e 7 -7) V e 7 / '"' 7(e 7 -7) 

e 7 W 7 2 (e 7 -7)^ 



1 \ e 4 
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where 7 = E e ±/m c c 2 , e 7 = i? 7 /m c c 2 , and e p h = E p h/m c c 2 . A numerical integration yields 
a pair spectrum that can be well fitted by a power law iV e ± oc E~2' 9 - The proportionality 
constant K 7y ^ e ± can be obtained from the absorbed gamma-ray luminosity. 

The presence of a magnetic field in the jet will imply that all these secondary pairs 
will produce synchrotron emission. Following Bosch-Ramon et al. (2005b) we assume that 
the magnetic field is entangled to cold protons in such a way it has random directions and 
hence the synchrotron emission is isotropic in the jet's frame. To calculate the synchrotron 
luminosity we estimate the specific emission (j e (z)) and absorption (k € (z)) coefficients from 
the secondary particle distribution (see Pacholczyk 1970 for the detailed formulae), in such 
a way that: 



= 2 ^W) X [1 " expHjA; ^ ))] ' (23) 

where to simplify the notation we are not using now primes to indicate that the calculation 
is in the jet's frame. In Eq. (23) Zj ~ i?j is the typical size of the synchrotron emitting 
plasma and e is the photon energy in units of m c c 2 . Integrating over the jet length we get 
the spectral energy distribution as: 

/*^-max 

obs 



T oas = f 

syn 

'20 



L (24) 



where 5 is the Doppler boosting factor defined as: 

6 = r(l-/3 b cos^ obs )' (25) 

To calculate the specific emission j e (z) we adopt different values of the magnetic field 
at z : Bq = 1, 10, and 100 Gauss (Bosch- Ramon & Paredes 2004). In Figure 5 we show 
the spectral energy distribution of the synchrotron radiation of all secondary pairs for the 3 
different values of Bq. The radio emission is quite negligible in comparison to the observed 
values, which at the minimum imply a luminosity of ~ 10 31 erg s _1 (e.g. Ribo et al. 2005). 



6. Discussion 

The predicted gamma-ray luminosity is clearly at its maximum during the periastron 
passage, when the neutron star travels through the densest parts of the wind. This is in 
accordance with the fact noticed by Massi (2004) that the peaks of the EGRET flux are 
coincident with the periastron and not with the radio maxima. The radio outbursts are 
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the result of particle injection in the jet that occurs after some relaxation time from the 
periastron passage, when the accretion rate is increased (Paredes et al. 1991). Any purely 
leptonic model for the gamma-ray emission would have to explain why the radio and gamma- 
ray peaks are not observed in similar orbital phases. 

Other specific feature of the gamma-ray emission predicted by our model is the presence 
of a local, secondary maximum at ip ~ 0.65 when the accretion rate, given by (1), has also 
a local maximum due to the fact that the wind velocity is roughly parallel to the neutron 
star orbital velocity, hence reducing v TC \ and increasing M c , as noticed by Marti & Paredes 
(1995). 

The effects of the opacity of the ambient photon fields to gamma-ray propagation pro- 
duces a "valley" in the spectral energy distribution, between a few tens of GeV and a few 
TeV, with a local minimum at around 100 GeV, during the periastron passage. The pre- 
dicted luminosity is within the detection possibilities of an instrument like MAGIC, which, 
integrating over several periastron passages, could build up a SED which can be compared 
with that presented in Fig. 3. Upper limits obtained with the Whipple telescope (Hall et 
al. 2003, Fegan et al. 2005) are indicated in the figure. The source is too weak for the 
sensitivity of this instrument according to our model. 

7. Concluding remarks 

We have presented a hadronic model for the high-energy gamma-ray production in the 
microquasar LS I +61 303. The model is based on the interaction of a mildly relativistic jet 
with a small content of relativistic hadrons with the dense equatorial disk of the companion 
B0 V star. Gamma-rays are the result of the decay of neutral pions produced by pp collisions. 
Charged pion decay will lead to neutrino production, that will be discussed elsewhere. The 
model takes into account the opacity of the ambient photon fields to the propagation of the 
gamma-rays. The predictions include a peak of gamma-ray flux in the periastron passage, 
with a secondary maximum at phase ip ~ 0.65. The spectral energy distribution presents a 
minimum around 100 GeV due to absorption. The spectral features should be detectable by 
an instrument like MAGIC through exposures ~ 50 hr, integrated along different periastron 
passages. 
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Table 1: Basic parameters assumed for the model 



Parameter 


Symbol 


Value 


Mass of the compact object 


M c 


1.4 M Q 


Jet's injection point 


z 


50 


Initial radius 


R 


zo/10 


Mass of the companion star 


M* 


10 M 


Radius of the companion star 


R* 


10 R Q 


Effective temperature of the star 


T e s 


22500 K 


Density of the wind at the base 


Po 


10~ n gr cm -3 


Initial wind velocity 


v 


5 km s _1 


Jet's Lorentz factor 


r 


1.25 


Minimum proton energy 


771 /min 

E p 


1 GeV 


Maximum proton energy 


7-1 /max 


100 TeV 


Penetration factor 


u 


0.1 


Eccentricity 


e 


0.72 


Orbital period 


P 


26.496 d 


Index of the jet proton distribution 


a 


2.2 


1 R g = GM c /c 2 . 
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Fig. 1. — Sketch of the situation described in the paper (not drawn to scale). 
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Fig. 2. — Orbit of LS I +61 303 drawn to scale illustrating the geometry used in the calcu- 
lations. The periastron passage occurs at 0.23 (Casares et al. 2005). 
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Fig. 3. — Upper, Left: A three-dimensional plot shows the generated luminosity as function 
of the orbital phase and gamma-ray energy. Right: The same plot, taken into account the 
gamma-ray absorption in the ambient photon fields. Down, Left: Ligthcurve for gamma- 
rays of energy 100 GeV. The dashed curve corresponds to the generated luminosity, whereas 
the continuous curve takes into account the effects of photospheric opacity. Right: Spectral 
energy distribution at the periastron and apastron passage. The unabsorbed spectra are in 
dashed lines. Upper limits from Whipple observations are indicated. 
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Fig. 4. — Opacity of the ambient photon fields to gamma-ray propagation. The figure 
corresponds to a viewing angle # b s = 30° and ip = 0.23 (periastron) . 
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